% 测试基矢产生

clearvars
sites = 13;
N_total = 13;

tic
memoization_matrix = cell(sites,N_total); 
[basis_matrix_new,~] = generate_basis_new(sites,N_total,memoization_matrix);
basis_matrix_new = sparse(basis_matrix_new); % 避免用==识别ID失败
clear('memoization_matrix')

fprintf('Dimension of the basis is %d\n',size(basis_matrix_new,1))
toc

% %% 对比和原版的区别
% memoization_matrix = cell(sites,N_total); 
% [basis_matrix,~] = generate_basis(sites,N_total,memoization_matrix);
% clear('memoization_matrix')
% 
% % 根据格点数产生质数序列，质数一定是无平方因子的数
% big_n = 1000;
% prim_numbers = primes(big_n);
% while length(prim_numbers)<sites
%     big_n = big_n*10;
%     prim_numbers = primes(big_n);
% end
% prim_numbers = prim_numbers(1:sites);
% % 根据质数对各个基矢产生唯一ID
% gene_ID_vector = sqrt(prim_numbers)';
% 
% IDs = basis_matrix*gene_ID_vector;
% IDs_new = basis_matrix_new*gene_ID_vector;
% logi = (IDs==IDs_new);
% find(logi~=1)
% % a=full(basis_matrix_new);
% % b=full(basis_matrix);
% % v1 = basis_matrix(23,:);
% % v2 = basis_matrix_new(23,:);
% % v3 = basis_matrix([23,24],:);
% % v4 = basis_matrix_new([23,24],:);
% % disp(v1==v2)
% % disp(v1*gene_ID_vector==v2*gene_ID_vector)
% % disp(v3==v4)
% % disp(v3*gene_ID_vector==v4*gene_ID_vector)
% % a(23,:)==b(23,:)
% % a(23,:)*gene_ID_vector==b(23,:)*gene_ID_vector
% % sum(abs(IDs-IDs_new))
%--------------------------------------------------------------------------------------------------
function [basis_matrix,memoization_matrix] = generate_basis_new(sites,N_total,memoization_matrix)
% fprintf('current problem is %d,%d\n',sites,N_total)
% 递归程序，用来产生占有数矢量（互相正交，但是这个矢量在MATLAB中直接内积不为0）
% memoization_matrix用来避免重复的计算，这个不需要应用于N_total==0和sites == 1的情况，因为结果已经有了
if N_total==0
    %如果没得放了就剩下的全补0
    % 全补充0的话用sparse会比用full矩阵快得多，nice，这样少粒子极限一万个格点一个粒子就可以瞬间产生了
    basis_matrix = sparse(zeros(1,sites));
    return
end
if sites == 1
    %如果最后一个site就剩下的全放上去
    basis_matrix = N_total; 
    return
end
% 有大量重复的问题：如果是已经解决的问题，记录下来，下次遇到就不要重新计算了
if isempty(memoization_matrix{sites,N_total})==0
    basis_matrix = memoization_matrix{sites,N_total};
    return
end

% 一般情况
% 下面都是为了提前分配好内存所做的
D=prod((N_total+1:N_total+sites-1)./(1:sites-1));
% 浮点计算非常接近这个整数，而不是别的整数，下面验证之，这个精度可以调
if abs(round(D)-D)>1e-8
    error('D not near integer')
end
D = round(D);
% 多格点少粒子极限用sparse才快
if sites<20
    basis_matrix = zeros(D,sites);
else
    basis_matrix = sparse(D,sites);
end
current_row = 1;
for n=0:N_total %选择当前第一个site填占据数n
    [sub_matrix,memoization_matrix] = generate_basis_new(sites-1,N_total-n,memoization_matrix);
    full_sub_matrix = [ones(size(sub_matrix,1),1)*n,sub_matrix]; 
    % 加列用逗号的
    sub_D = size(full_sub_matrix,1);
    basis_matrix(current_row:current_row+sub_D-1,:) = full_sub_matrix;
    current_row = current_row+sub_D;
end
memoization_matrix{sites,N_total}=basis_matrix;

% 如果memoization_matrix太多，及时清空以前的，以免耗尽内存
% 清理的规则来自对运行的观察，格点数会从小变大。放到basis_matrix赋值之后不影响正常运行
if rem(sites,10)==0
    for i=sites-10+1:sites-5+1
        for j=1:size(memoization_matrix,2)
            memoization_matrix{i,j}=[];
        end
    end
end
end
%% new是最终版本，下面的版本是用来对比是否正确的
%% 太好了，我自己的算法成功超越了European Journal of Physics, 2010, 31(3): 591.中的算法

%-------------------------------------------------------------------------------------------

% 下面是不提前分配内存的函数，用于和上面的对比
function [basis_matrix,memoization_matrix] = generate_basis(sites,N_total,memoization_matrix)
% fprintf('current problem is %d,%d\n',sites,N_total)
% 递归程序，用来产生占有数矢量（互相正交，但是这个矢量在MATLAB中直接内积不为0）
% memoization_matrix用来避免重复的计算，这个不需要应用于N_total==0和sites == 1的情况，因为结果已经有了
if N_total==0
    %如果没得放了就剩下的全补0
    % 全补充0的话用sparse会比用full矩阵快得多，nice，这样少粒子极限一万个格点一个粒子就可以瞬间产生了
    basis_matrix = sparse(zeros(1,sites));
    return
end
if sites == 1
    %如果最后一个site就剩下的全放上去
    basis_matrix = N_total; 
    return
end
% 有大量重复的问题：如果是已经解决的问题，记录下来，下次遇到就不要重新计算了
if isempty(memoization_matrix{sites,N_total})==0
    basis_matrix = memoization_matrix{sites,N_total};
    return
end

% 一般情况
basis_matrix = [];
for n=0:N_total %选择当前第一个site填占据数n
    [sub_matrix,memoization_matrix] = generate_basis(sites-1,N_total-n,memoization_matrix);
    full_sub_matrix = [ones(size(sub_matrix,1),1)*n,sub_matrix]; 
    % 加列用逗号的
    basis_matrix = [basis_matrix;full_sub_matrix];
end
memoization_matrix{sites,N_total}=basis_matrix;

% 如果memoization_matrix太多，及时清空以前的，以免耗尽内存
% 清理的规则来自对运行的观察，格点数会从小变大。放到basis_matrix赋值之后不影响正常运行
if rem(sites,20)==0
    for i=sites-20+1:sites-10+1
        for j=1:size(memoization_matrix,2)
            memoization_matrix{i,j}=[];
        end
    end
end

end